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The  coving  finite  element  (MFE)  method  is  a  new  PDF  solution  method  which  has 
shown  significant  promise  in  1-D  for  the  numerical  solution  of  some  of  the  mos 
difficult  problems  under  study  with  extremely  large,  but  finite,  gradients. 

The  overall  objective  of  the  present  research  is  to  explore  further  the  promis 
of  the  continuous  node  moving  properties  of  the  MFE  method  in  2-D.  For  this, 
both  the  logical  structure  of  the  KFE  method  and  its  reduction  to  practice  in 
2-D  are  under  investigation  in  this  project.  This  initial  CONTINUED ) 
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I.  MATHEMATICAL  INTRODUCTION 


In  order  to  discuss  concisely  the  progress  of  research  during  this 
reporting  period,  the  basic  formulation  of  the  MFE  method  in  2-0  is 
presented  immediately  below. 

Consider  a  general  system  of  PDE's,  U  =  L(U),  or 

ui  =  L ! ( U ) 


UN  =  ln(u) 


(1) 


Using  piecewise  linear  approximants  of  uk  ~  ^  akj^  •  a*<J(x,  y) 

J 

on  a  hexagonally  connected  triangular  mesh  (see  Figure  1),  application  of 
the  chain  rule  to  the  differentiation  of  u^  gives 


u.  =  >:  ak.ukJ  +  x.nkJ  +  y.ykJ 
K  i  J  J  J 


,  where 


i  uUk 
ukJ  = 

uakj 


j  i)uk  i  3uk 
pkJ  =  ---  ,  >kJ  =  ; — — 

<)Xj  Jyj 


(2) 


(3) 


The  functions  ak^ ,  gk1 ,  and  uk3  are  basis  functions  at  each  node,  j, 
corresponding  to  the  effect  on  the  piecewise  linear  functions  uk  by  the  N 
+  2  parameters  akj,  xj,  and  yj.  It  can  be  shown  readily  by  taking 
tiny  increments  in  ak j ,  or  xj,  or  yj  that  all  three  basis  functions 
«,  a  and  y  are  piecewise  linear  functions  having  their  support  in  the 
hexagon  of  six  triangles  surrounding  the  jth  node.  It  turns  out  that 
<«kJ  is  independent  of  k;  it  is  continuous  and  takes  on  the  value  1  at  the 
center  of  the  hexagon  and  0  at  the  other  nodes.  The  functions  and 
are  discontinuous  at  the  center  and  on  the  inner  edges  of  the  hexagon;  they 
vanish  identically  on  the  outer  edges  of  the  hexagon. 

Ordinary  differential  equations  (ODE's)  are  derived  for  the  parameters 
of  the  MFE  method  by  requiring  that  the  parameter  derivatives  {al..,  ...aN.,  x.,  y.} 
are  evaluated  at  each  instant  so  as  to  formally  minimize  the  L2  norm  of 
the  POE  residuals,  U  -  L(U),  plus  regularization  terms.  The  variational 
equations  for  this  minimization  yield  the  following  system  of  ODE's: 


-1 


u 


FIGURE  1-A.  Exact  solution  surface,  with  lines  of 

CONSTANT  X  AND  CONSTANT  Y, 


APPROXIMATION  SURFACE 

U*  55  F*(X,Y) 


FIGURE  1-B.  Approximate  solution  represented  by 

PIECEWISE  LINEAR  FUNCTIONS  MAKING  UP 
HEXAGONALLY  CONNECTED  TRIANGULAR  FACETS. 
NFE  UNKNOWNS  ARE  lal j  .  •  • . -aNj  ,  x  . ,  y  ■ ) 


J)ak.  +  (tikJ ,  a1);.  +  (ykJ,  a1  )y  •  =  (L  (U),  a1)  for  k  =  1,  ...  N,  (4a) 
J  J  J  * 

<«j.  iik'i)akj  +  iik’)ij  +  (yk^,  pk’  )y^ 

N  i 

=  >;  (L. (U),  3k  )  +  (regularization  terms)  (4b) 

k=l 

;uj,  Yk^Jakj  +  (gkj,  yk’jxj  +  (yk^,  yk1)^ 

N  i 

=  z  (L.  (U),  yk  )  +  (regularization  terms)  (4c) 
k=l  K 


j 

j 


The  sums  on  j  in  Eqns.  (4)  run  over  the  seven  neighboring  nodes  of 
i  (including  the  i^h  node  itself)in  the  hexagonal  grid.  Equations  (4)  thus 
provide  the  basic  working  equations  of  the  MFE  method  in  2-D.  This  system 
of  ODE'S  is  written  for  purposes  of  numerical  solution  in  the  form 


C(y)y  -  g(y)  =  0  ,  (5) 

where  y(t)  =  (al^  ....  yj;  al2,  ....  x2,  y2l  . )  is 


the  vector  of  unknown  parameters,  and  the  "mass  matrix"  C(y)  is  symmetric  j 

and  positive  definite.  This  system  of  ODE's  can  be  extremely  stiff,  and  j 

such  stiff  ODE  solvers  as  the  Gear  method  and  an  implicit  Runge-Kutta  method  < 

are  used  for  the  numerical  solutions  of  Equations  (4)  and  (5).  ] 

Obviously,  the  evaluation  of  the  numerous  inner  product  terms  in  Eqns. 

(4)  is  tedious,  and  the  overall  organization  of  an  efficient  code  structure 
for  the  general  solution  of  Eqns.  (4)  and  (5)  in  2-D  is  quite  exacting. 

These  are  the  major  tasks  which  have  been  performed  during  this  reporting 
period  and  which  will  be  summarized  and  discussed  further  below. 

J 

\ 


(  j 

•  Vu  , 

j 

N 

j:  z 

k=l  j 


N 

x  j: 

k=l  j 
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RESEARCH  PROGRESS 


The  overall  otr  'tive  of  this  research  is  to  explore  the  promise  of  the 
continuous  node  moving  properties  in  2-0  of  the  MFE  method  (a  new  numerical 
POE  solution  method)  for  the  effective  solution  of  problems  with  extremely 
large,  but  finite,  gradients.  Such  problems  arise  in  numerous  scientific 
and  engineering  applications  which  are  of  interest  to  the  U.S.  Air  Force. 
Accordingly,  this  research  necessarily  involves  both  fundamental  analysis  of 
the  njnerical  method,  itself,  and  reduction  to  practice  of  the  MFE  method 
in  order  to  carry  out  the  essential  scientific  testing  and  advancement  of 
the  basic  method  in  higher  dimensions. 

Problems  which  can  be  written  in  conservation  form  as 

ut  =  -Ci.fx  ~  c2‘9y  +  c3’(uxx  +  uyy)  (6) 

are  being  investigated  in  the  current  stages  of  research.  The  immediate 
objective  of  this  first  year’s  research  is  to  develop  a  working  MFE  code  in 
2-D  which  can  solve  such  simple  equations  as  the  heat  equation,  a  single 
travelling  wave,  and  Burger's  equations,  in  order  to  study  at  an  early  date 
the  basic  node  moving  properties  of  the  MFE  method  in  2-D.  (All  of  these 
simple  equations  are  obtained  by  appropriate  choices  of  coefficients  cj, 

C2»  and  C3  and  functions  fx  and  gy  in  Eqn.  (6).)  The  first  working 
code  version  is  nearly  operational  at  this  time.  Table  1  presents  a  brief 
summary  of  the  status  of  essential  tasks  conducted  in  this  research. 


DISCUSSION  AND  ADDITIONAL  INFORMATION 


The  work  during  this  reporting  period  has  been  heavily  devoted  to  code 
development  and  analysis  of  the  logical  structure  of  the  MFE  method  in  two 
dimensions.  As  a  result  of  these  efforts,  significant  computational 
economies  have  been  realised  for  both  unvectorized  versions  of  the  method 
as  it  currently  exists  and,  more  importantly,  for  vectorized  versions  which 
will  undoubtedly  emerge  in  later  efforts. 
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TABLE  1  -  SUMMARY  OF  PROGRESS 


Basis  Functions 


MILESTONE:  Using  piecewise  linear 
basis  functions,  develop  inner 
products  of  the  operators  in  Eqn. 
(6)  by  nunerical  integration. 


Regularization 

MILESTONE:  Develop  penalty  functions 
which  prevent  grid  mesh  tangling  for 
the  MFE  nodes. 

Grid  Mesh  Generation 

MILESTONE:  Assess  automatic  grid 
mesh  generators  for  compatibility 
with  the  data  structures  in  the 
MFE  method. 

Matrix  Solution 

MILESTONES: 

(a)  Apply  existing  direct  factori¬ 
zation  methods. 

(b)  L-U  decomposition  of  banded 
matrix. 

(c)  Point  and  line  relaxation  methods. 

(d)  Alternating  Direction  Implicit 
(ADI). 


Boundary  Conditions 
MILESTONES: 

(a)  Stationary  Dirichlet  and  zero- 
Neumann  conditions. 

(b)  Time-dependent  Dirichlet 
conditions. 

(c)  Non-zero  Neumann  and  mixed 
boundary  conditions. 


STATUS:  Completed.  Some  operators 
in  gas  dynamics  equations  have  also 
been  evaluated,  coded,  and  tested. 
Both  simple  and  composite  Simpson's 
(Newton-Cotes)  methods  are  used  on 
triangle  edge  (line)  integrations. 
Both  simple  and  composite  midpoint 
rules  are  used  for  integrations  on 
triangle  areas. 


STATUS:  Completed. 


STATUS:  Scoping  completed. 


STATUS:  Completed. 

STATUS:  Completed. 

STATUS:  Completed. 

STATUS:  Two-directional  ADI 
completed.  Dynamic  ADI  in  three 
directions  warrants  additional 
research,  as  do  efficient  iterative 
matrix  solution  methods. 


STATUS:  Completed. 
STATUS:  Completed. 

:  Barely  started. 


STATUS 


The  nature  of  these  economies  are  illustrated,  for  example,  by  the 
structure  of  equations  (4)  above  which  are  written  on  a  node-by-node  basis 
(via  the  s  unimat  ions  over  j  on  nodes).  Although  compact  for  purposes  of 
theoretical  derivations  and  expositions,  this  node-by-node  ordering  of 
operations  turns  out  not  to  provide  the  most  efficient  architecture  for  MF E 
computations  of  those  PDE  operators  which  appear  in  Eqn.  (6).  Alternatively, 
the  structure  of  the  MFE  method  and  the  topology  of  hexagonal ly  connected 
triangular  basis  functions  lend  themselves  to  both  logical  and  practical 
computational  economies  when  the  method  is  executed  on  a  triangle-by- 
triangle  basis.  This  alternative  course  is  the  one  which  is  being 
implemented  in  the  present  work. 

In  scoping  studies  of  grid  mesh  generation  methods,  it  has  become 
evident  that  hexagonal ly  connected  triangular  mesh  generation  methods  are 
basically  compatible  with  the  MFE  method.  The  logical  structure  and  mode  of 
PDE  solution  execution  by  the  MFE  method  again  contributes  to  this  conclu¬ 
sion.  So  long  as  the  MFE  method  continues  to  exhibit  robustness  in  its 
continuous  node  movement  properties,  the  use  of  grid  mesh  generation 
routines  would  be  required  only  for  problem  initializations.  All  solution 
regridding  subsequent  to  initialization  is  continuously  performed,  in 
effect,  by  the  MFE  method  itself  based  upon  its  own  stringent  minimization 
criteria.  It  is,  in  fact,  preferable  to  not  introduce  a  generically 
different  mesh  generation  method  on  an  in-line  basis  into  the  transient  MFE 
solutions.  This  anticipated  restriction  to  problem  initializations  reduces 
greatly  the  role  and  the  demands  upon  other  mesh  generation  methods  which 
would  be  used  in  conjunction  with  the  MFE  method.  In  those  problems  which 
may  require  some  infrequent  remapping  of  MFE  solutions  in  order  to  reach 
final  solutions,  one  would  undoubtedly  choose  to  interrupt  the  MFE  solution, 
regrid  the  numerical  PDE  solution  data  at  that  stage,  and  then  proceed  again 
with  the  MFE  solution  as  a  new  initial-value  problem.  For  this  type  of 
procedure,  the  MFE  method  would  be  viewed  as  the  dynamic  mesh  generation 
procedure  and  the  alternative  mesh  generation  method  would  be  viewed  as  the 
static  in i t i a  1  izat  ion  procedure.  It  thus  appears  that  if  any  revisions  at 
all  are  required,  it  would  be  only  minor  alterations  of  mesh  indexing  and 
of  triangle  orientations  which  may  be  needed  in  some  of  the  available 
triangular  mesh  generation  packages  in  order  to  serve  effectively  as  an 
initialization  method  for  MFE  solutions.  Further  work  in  this  area  of  mesh 
generation  for  initialization  will  be  taken  up  again  in  later  stages  of 
research  when  more  complex  applications  are  considered. 


Matrix  solution  methods  for  the  large,  handed,  stiff  system  of  MFE 
equations  have  been  tested  extensively  during  this  reporting  period.  The 
results  of  this  testing  indicate  that  L  -  U  decomposition  of  the  banded 
matrix  of  the  MFE  method  is  a  reasonable  choice  for  this  early  stage  of 
research  where  immediacy,  reliability,  and  ease  of  implement  at  ion  outweigh 
computational  efficiency  in  small  test  problem  applications.  (In  later 
stages  of  research  it  will  be  necessary  to  implement  more  efficient  matrix 
solution  methods  for  larger  and/or  more  complex  POE  systems,  as  will  be 
indicated  in  the  discussion  below.)  Point  and  line  relaxation  methods  have 
also  been  tested  during  this  reporting  period  and  found  to  converge  too 
slowly  to  be  useful  in  current  MFE  research.  Two-d irect ional  ADI  methods 
were  similarly  tested  and  were  found  to  be  well-suited  for  implementation 
with  the  MFE  method  on  a  quadr i lateral  grid.  For  the  hexagonally  connected 
triangular  mesh  which  is  used  in  the  present  effort,  three-d irect ional  (in 
contrast  to  two-direct ional )  ADI  methods  will  be  required.  Although  a  large 
amount  of  additional  research  is  needed  before  three-direct ional  ADI  can  be 
used  effectively  with  the  MFE  method,  dynamic  ADI  methods  are  sufficiently 
promising  that  such  further  research  certainly  appears  to  be  warranted  for 
large-scale  MFE  applications  in  the  future .  Indeed,  dynamic  ADI  and  perhaps 
other  splitting  methods  presently  appear  to  be  mandatory  in  order  to  achieve 
desired  economies  when  large-scale  applications  are  considered.  In  antici¬ 
pation  of  these  later  matrix  solution  needs,  the  investigation  of  more 
efficient  interative  solution  methods  has  also  been  started.  At  this  time, 
partial  block  Gauss  elimination  and  incomplete  orthogonal i zat ion  of  a 
conjugate  gradient  solution  method  are  in  the  early  scoping  stages  and  are 
showing  some  signs  of  promise. 

No  major  deviations  from  the  planned  course  of  research  have  been  made 
during  this  reporting  period  and  none  are  anticipated  during  the  next 
reporting  period.  So  far  as  maintaining  an  alert  for  potential  problems 
which  may  arise,  the  adequacy  of  numerical  Jacobians  is  under  constant 
survei 1 1 ance.  Testing  during  the  present  period  indicates  that  essential 
levels  of  computational  accuracy  in  Jacobian  evaluations  are  likely  to  be 
maintained  by  the  numerical  discret izat ion  methods  which  are  used  in  current 
code  versions.  Constant  monitoring  of  this  portion  of  MFE  calculations  is 
warranted  because  such  a  large  and  tedious  coding  effort  of  analytic 
Jacobian  expressions  would  be  required  in  the  event  that  numerical  Jacobian 
evaluations  were  found  to  be  inadequate. 
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